//Burgers--RK2
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.1415926

extern inline double ux(double x);
extern inline double fux(double ux);

typedef struct geo_mesh
{
    double u;        //state value of the next step
    double u0;       //state value of the present step
    double xc;       //cell center coordinate
    double res_buff; //only for RK3
    double res;      //f[i+1/2]-f[i-1/2]
    double du;       //only for reconstruction
    double k;        //only for reconstruction
} geo;

int ncells;
geo *cell; //pointer for memory
double dx;

/**
 * @brief get dt for each time step. Named Courant-Friedrichs-Lewy condition, 
 * CFL=c*dt/dx, here c represents the speed of the wave.
 * @param cfl 
 * @param dx 
 * @param u(x)
 * @return dt
 */
double timestep(double cfl, double dx, int ncells)
{
    double c = -1.0; //let c = -1
    double temp;
    for (int i = 1; i <= ncells; i++)
    {
        temp = fabs(cell[i].u);
        if (c > temp)
            c = c;
        else
            c = temp;
    }
    return cfl * dx / c;
}

extern inline double ux(double x)
{
    return (1 - cos(2 * pi * x));
}

extern inline double fux(double x)
{
    return 0.5 * x * x;
}

/**
 * @brief limiter
 * @param k_left
 * @param k_right
 * @return min_k
 */
double minmod(double a, double b)
{
    double k;
    if (fabs(a) < fabs(b))
        k = a;
    else
        k = b;
    return k;
}

/**
 * @brief upwind 2 with limiter minmod
 * @param cell[i].u0
 * @param cell[i].dul
 * @param cell[i].dur
 * @param cell[i].du
 * @return interface flux for res 
 */
double flux(double u, double v, double w)
{
    double flx;
    double dul;
    double dur;
    double du;
    double k;
    dul = fabs(v - u);
    dur = fabs(w - v);
    du = minmod(dul, dur); //get du
    k = du / dx;
    flx = fux(v + dx * 0.5 * k);
    return flx;
}

/**
 * @brief output file in Tecplot .plt
 * @param ncells 
 * @param N No. of the file
*/
void output(int ncells, int N, const char *zonename, int i, int ii)
{
    FILE *fp;
    char fname[64];
    if (i % ii == 0)
    {
        sprintf(fname, "solution_%04d.plt", N);
        fp = fopen(fname, "w");
        fprintf(fp, "variables=\"x\",\"u\"\n");
        fprintf(fp, "Zone T=\"%s step---%05d\"\n", zonename, N);
        for (int i = 1; i <= ncells; i++)
        {
            fprintf(fp, "%lf %lf \n", cell[i].xc, cell[i].u); //output with the center of cell
        }
        fclose(fp);
    }
    return;
}

/**
 * @brief f[i+1/2]-f[i-1/2]
 * @param flux
 * @return res
 */
void res_cal(void)
{
    int i;
    for (i = 1; i <= ncells; i++)
        cell[i].res = 0; //init per step
    for (i = 2; i <= ncells - 1; i++)
        //inner cell
        cell[i].res = flux(cell[i - 1].u, cell[i].u, cell[i + 1].u) - flux(cell[i - 2].u, cell[i - 1].u, cell[i].u);
    //boundary cell
    cell[1].res = flux(0.0, cell[1].u, cell[2].u);
    cell[ncells].res = flux(cell[ncells - 1].u, cell[ncells].u, 0.0);

    return;
}

/**
 * @brief time advancing 
 * @param res
 * @return next step u(x)
 */
void rk2(int ncells, double dt, double dx)
{
    for (int istage = 1; istage <= 2; istage++) // 2 stages
    {
        res_cal();       //get f[i+1/2]-f[i-1/2] for step 1 and step 2 respectivly
        if (istage == 1) //solution update
            for (int j = 1; j <= ncells; j++)
            {
                cell[j].u0 = cell[j].u;                          //save current results
                cell[j].u = cell[j].u - (dt / dx) * cell[j].res; //get u*
            }
        else
            for (int j = 1; j <= ncells; j++)
            {
                cell[j].u = cell[j].u - (dt / dx) * cell[j].res;
                cell[j].u = 0.5 * (cell[j].u0 + cell[j].u); //get next step u
            }
    }
    return;
}

int main()
{
    char const *case_name = "Burgers_RK2";
    double xmin, xmax, t;
    double tf = 5.0; //final time
    double cfl = 0.8;
    double dt;
    int cfl_condition;
    int i, ii;
    int nsteps;
    int itime; //control max time steps
    xmin = 0;
    xmax = 10;
    ncells = 2000;
    dx = (xmax - xmin) / ncells;
    cell = calloc(ncells + 2, sizeof(geo)); //2 ghost cells

    printf("time step switch, press :\n \"0 \"for fixed dt = 0.002\n \"1\" for cfl\n");
    scanf("%d", &cfl_condition);
    printf("Output per step:\n");
    scanf("%d", &ii);

    if (cfl_condition)
        printf("CFL condition applied, for each step:\n");
    //initial condition
    for (i = 0; i <= ncells + 1; i++)
    {
        cell[i].xc = 0.025 + (i - 1) * dx; //get the ceter coordinate of the cell
        if (cell[i].xc >= 0.0 && cell[i].xc <= 0.975)
            cell[i].u0 = ux(cell[i].xc); //replace u0 with average data
        else
            cell[i].u0 = 0.0;
        cell[i].u = cell[i].u0; //1 step
    }

    //time advance
    t = 0.0;
    nsteps = 0;
    for (itime = 0; itime <= 5000; itime++) //general control, max step <= 5000
    {
        //output current solution
        output(ncells, itime, case_name, itime, ii);
        if (t == tf)
            break;
        //choose time step
        if (cfl_condition)
        {
            double dt_rk2 = timestep(cfl, dx, ncells);
            dt = dt_rk2;
        }
        else
            dt = 0.002;
        if (t + dt > tf)
            dt = tf - t;
        t = t + dt;
        nsteps = nsteps + 1;
        if (itime % ii == 0)
            printf("dt=%lf, step=%d\n", dt, nsteps - 1);
        rk2(ncells, dt, dx);

        for (i = 0; i <= ncells + 1; i++)
            cell[i].u0 = cell[i].u; //update
    }
    //output(ncells, -1, case_name);
    printf("\nfinal time t(sec) = %20.3f by %5d time steps\n", t, nsteps);

    free(cell); //deallocate cell
    cell = NULL;

    return 0;
}
